home *** CD-ROM | disk | FTP | other *** search
/ The Atari Compendium / The Atari Compendium (Toad Computers) (1994).iso / files / prgtools / mint / lib / mntlib44.zoo / mntlib / _mulsf3.cpp < prev    next >
Encoding:
Text File  |  1993-06-17  |  5.1 KB  |  225 lines
  1. | mjr: not needed on the TT
  2.  
  3. #ifndef    __M68881__
  4.  
  5.     .text
  6.     .even
  7.     .globl    __mulsf3, ___mulsf3
  8.  
  9. __mulsf3:
  10. ___mulsf3:
  11.  
  12. # ifdef    sfp004
  13.  
  14. | single precision floating point stuff for Atari-gcc using the SFP004
  15. | developed with gas
  16. |
  17. | single floating point multiplication routine
  18. |
  19. | M. Ritzert (mjr at dfg.dbp.de)
  20. |
  21. | 7. Juli 1989
  22. |
  23. | revision 1.1: adapted to the gcc lib patchlevel 58
  24. | 4.10.90
  25.  
  26. comm =     -6
  27. resp =    -16
  28. zahl =      0
  29.  
  30.     lea    0xfffffa50:w,a0
  31.     movew    #0x4400,a0@(comm)    | load first argument to fp0
  32.     cmpiw    #0x8900,a0@(resp)    | check
  33.     movel    sp@(4),a0@
  34.     movew    #0x4427,a0@(comm)
  35.     .long    0x0c688900, 0xfff067f8
  36.     movel    sp@(8),a0@
  37.     movew    #0x6400,a0@(comm)    | result to d0
  38.     .long    0x0c688900, 0xfff067f8
  39.     movel    a0@,d0
  40.     rts
  41.  
  42. # else    sfp004
  43.  
  44. | single floating point multiplication routine
  45. |
  46. | written by Kai-Uwe Bloem (I5110401@dbstu1.bitnet).
  47. | Based on a 80x86 floating point packet from comp.os.minix, written by P.Housel
  48. |
  49. | Revision 1.2.4 michal 05-93 (ntomczak@vm.ucs.ualberta.ca)
  50. |   + ensure that Inf * NaN == NaN * Inf == NaN
  51. |     and 0 * Inf = Inf * 0 = NaN
  52. |
  53. | Revision 1.2.3 michal 05-93 (ntomczak@vm.ucs.ualberta.ca)
  54. |   + code smoothing
  55. |
  56. | patched by Olaf Flebbe (flebbe@tat.physik.uni-tuebingen.de)
  57. |
  58. | Revision 1.2.2 olaf 05-93
  59. |   + fixed a bug with -0.
  60. |
  61. | Revision 1.2.1 olaf 12-92:
  62. |   + added support for NaN and Infinites
  63. |   + added support for -0
  64. |
  65. | Revision 1.2, kub 01-90 :
  66. | added support for denormalized numbers
  67. |
  68. | Revision 1.1, kub 12-89 :
  69. | Created single float version for 68000. Code could be speed up by having
  70. | the accumulator in the 68000 register set ...
  71. |
  72. | Revision 1.0:
  73. | original 8088 code from P.S.Housel for double floats
  74.  
  75. BIAS4    =    0x7F-1
  76.  
  77.     lea    sp@(4),a0
  78.     moveml    d2-d5,sp@-
  79.     moveml    a0@,d4/d5    | d4 = v, d5 = u
  80.  
  81.     movel    #0x7fffff,d3
  82.     movel    d5,d0        | d0 = u.exp
  83.     andl    d3,d5        | remove exponent from u.mantissa
  84.     swap    d0
  85.     movew    d0,d2        | d2 = u.sign
  86.  
  87.     movel    d4,d1        | d1 = v.exp
  88.     andl    d3,d4        | remove exponent from v.mantissa
  89.     swap    d1
  90.     eorw    d1,d2        | d2 = u.sign ^ v.sign (in bit 15)
  91.     
  92.     moveq    #15,d3
  93.     bclr    d3,d0        | kill sign bit
  94.     bclr    d3,d1        | kill sign bit
  95.     tstl    d0        | test if one of factors is 0
  96.     beq    1f
  97.     tstl    d1
  98. 1:    seq    d2        | 'one of factors is 0' flag in the lowest byte
  99.     lsrw    #7,d0        | keep here exponents only
  100.     lsrw    #7,d1
  101.  
  102. |
  103. | Testing for NaN and Infinities
  104. |
  105.     moveq    #-1,d3
  106.     cmpb    d3,d0
  107.     beq    0f
  108.     cmpb    d3,d1
  109.     bne    nospec
  110.     bra    1f
  111. |
  112. |    first operand is special
  113. |
  114. 0:    tstl    d5        | is it NaN?
  115.     bne    retnan
  116. 1:    tstb    d2        | 0 times special or special times 0 ?
  117.     bne    retnan        | yes -> NaN
  118.     cmpb    d3,d1        | is the other special ?
  119.     beq    2f        | maybe it is NaN
  120. |
  121. |    Return Infinity with correct sign
  122. |
  123. retinf: movel    #0xff000000,d0  | we will return #0xff800000 or #0x7f800000
  124.     lslw    #1,d2
  125.     roxrl   #1,d0        | shift in high bit as given by d2
  126. return:    moveml    sp@+,d2-d5
  127.     rts
  128.  
  129. |
  130. | v is special
  131. |
  132.     
  133. 2:    tstl    d4        | is this NaN?
  134.     beq    retinf        | we know that the other is not zero
  135. retnan: moveql    #-1,d0
  136.     lsrl    #1,d0        | 0x7fffffff -> d0
  137.     bra    return
  138. |
  139. | end of NaN and Inf.
  140. |
  141. nospec:    tstb    d2        | not needed - but we can waste two instr.
  142.     bne    retzz        | return signed 0 if one of factors is 0
  143.     moveq    #23,d3
  144.     bset    d3,d5        | restore implied leading "1"
  145.     subqw    #8,sp        | multiplication accumulator
  146.     tstw    d0        | check for zero exponent - no leading "1"
  147.     bne    1f
  148.     bclr    d3,d5        | remove it
  149.     addqw    #1,d0        | "normalize" exponent
  150. 1:    tstl    d5
  151.     beq    retz        | multiplying zero
  152.  
  153.     moveq    #23,d3
  154.     bset    d3,d4        | restore implied leading "1"
  155.     tstw    d1        | check for zero exponent - no leading "1"
  156.     bne    1f
  157.     bclr    d3,d4        | remove it
  158.     addqw    #1,d1        | "normalize" exponent
  159. 1:    tstl    d4
  160.     beq    retz        | multiply by zero
  161.  
  162.     addw    d1,d0        | add exponents,
  163.     subw    #BIAS4+16-8,d0    | remove excess bias, acnt for repositioning
  164.  
  165.     clrl    sp@        | initialize 64-bit product to zero
  166.     clrl    sp@(4)
  167.  
  168. | see Knuth, Seminumerical Algorithms, section 4.3. algorithm M
  169.  
  170.     movew    d4,d3
  171.     mulu    d5,d3        | mulitply with bigit from multiplier
  172.     movel    d3,sp@(4)    | store into result
  173.  
  174.     movel    d4,d3
  175.     swap    d3
  176.     mulu    d5,d3
  177.     addl    d3,sp@(2)    | add to result
  178.  
  179.     swap    d5        | [TOP 8 BITS SHOULD BE ZERO !]
  180.  
  181.     movew    d4,d3
  182.     mulu    d5,d3        | mulitply with bigit from multiplier
  183.     addl    d3,sp@(2)    | store into result (no carry can occur here)
  184.  
  185.     movel    d4,d3
  186.     swap    d3
  187.     mulu    d5,d3
  188.     addl    d3,sp@        | add to result
  189.                 | [TOP 16 BITS SHOULD BE ZERO !]
  190.     moveml    sp@(2),d4-d5    | get the 48 valid mantissa bits
  191.     clrw    d5        | (pad to 64)
  192.  
  193.     movel    #0x0000ffff,d3
  194. 2:
  195.     cmpl    d3,d4        | multiply (shift) until
  196.     bhi    3f        |  1 in upper 16 result bits
  197.     cmpw    #9,d0        | give up for denormalized numbers
  198.     ble    3f
  199.     swap    d4        | (we''re getting here only when multiplying
  200.     swap    d5        |  with a denormalized number; there''s an
  201.     movew    d5,d4        |  eventual loss of 4 bits in the rounding
  202.     clrw    d5        |  byte -- what a pity 8-)
  203.     subqw    #8,d0        | decrement exponent
  204.     subqw    #8,d0
  205.     bra    2b
  206. 3:
  207.     movel    d5,d1        | get rounding bits
  208.     roll    #8,d1
  209.     movel    d1,d3        | see if sticky bit should be set
  210.     andl    #0xffffff00,d3
  211.     beq    4f
  212.     orb    #1,d1        | set "sticky bit" if any low-order set
  213. 4:    addqw    #8,sp        | remove accumulator from stack
  214.     jmp    norm_sf        | (result in d4)
  215.  
  216. retz:    addqw    #8,sp        | release accumulator space
  217. retzz:    moveq    #0,d0        | save zero as result
  218.     lslw    #1,d2        | and set it sign as for d2
  219.     roxrl    #1,d0
  220.     moveml    sp@+,d2-d5
  221.     rts            | no normalizing neccessary
  222.  
  223. # endif    sfp004
  224. #endif    __M68881__
  225.